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ABSTRACT 


The distribution of a measure of distance between a probability 
density function and its estimate is examined through Monte Carlo 
methods. The estimate of the density function is that proposed by 
Rosenblatt using sums of weight functions centered at the observed 
values of the random variables. The weight function used in all 
cases was triangular, but both uniform and Cauchy densities were 
tried for different sample sizes and bandwidths. The simulated dis- 
tributions appear in all cases to be close to Gamma distributions, 
but it has not been possible to relate the parameters to the popu- 


lation characteristics or to the window shape and bandwidth. 
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I. INTRODUCTION 


There are several recently proposed classes of empirical proba- 
bility density function [Refs. 1-5] all generally considered to be 
superior to the classical histogram estimates. The class considered 
in this paper is based on independent observations, i.e., Xie Xoreces 
x are independent and identically distributed random variables with 
continuous unknown density function f(x). The method used to esti- 
mate f(x) is that proposed by Rosenblatt [4] and is as follows. 


Denote the estimator of f(x) by £ (x) - Then 


1 n xh= X:. 
2S SET) zt b(n) } : 


~ 


where W(u) is a bounded integrable weight function with 


[> =] 
Sf W(u) du = 1 
Oo 
and b(n) is a bandwidth that tends to zero as n*>+©™,, but is such 


that o(b(n)) = aoe i.e. it converges at a rate slower than l1/n. 


Thus we might have b(n) ~ WAS 


- Note that all estimates of this 
form for given observations are themselves density functions; that 


is 
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Since the X.‘s are random variables, £ (x) is a continuous parameter 
stochastic process. 

The local properties of such estimates satisfying relatively madd 
conditions can be shown to be biased [3]. Our object in this thesis 
is to Wigest tigate o global measure of how good £ (x) is as an esti- 
mate of f(x). This measure was proposed by Bickel and Rosenblatt [2]. 
The particular global measure gisbieeel for investigation is the 


’ statistic: 


[E. (x) - £(x)1° 
B(n) - {= dx , 


whose value will vary with each realization of X coe, x, peeled 


1’ 
is a statistic or function of the n random variables. 
Since an exact distribution of this statistic is not mathematically 
tractable, we have examined some representative cases by simulation. 
The conjecture made by Bickel and Rosenblatt [2] is that the 
statistic B(n) is asymptotically normally distributed regardless of 
the underlying distribution f(x). With this conjecture in mind, we 
hoped to establish a criterion that can be used to perform goodness 
of fit testing. 
A problem with this type of criterion is that, unlike the Kolmo- 
Grove sei choy or Cramer-von-Mises statistics for goodness of fit of 
an empirical distribution function to the unknown distribution, B(n) 
is not distribution free and its distribution for finite n is unknown 
for any case. However, it is believed that B(n) is fairly robust 
with rapid convergence to the asymptotic distribution; another object 


of the thesis is to examine this hypothesis. 
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Going beyond asymptotic normality of B(n), the following result 
for a quadratic functional is also of some interest [2]. The func- 
tion a(x) is assumed to be a bounded piece-wise smooth integrable 
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{b(n) } tf tnb(n) f te, (x)~ £(x)]) a(x)dx - J £(x)a(x)dxf W(z) dz] 
is asymptotically normally distributed with mean zero and variance 


aw'®) (0) few? £ (x)? dx 


(4) 


as n*+>© , where W is the 4-th convolution of W with itself. 


ve galwotted eas stab 36 Xb ri fata 

anu) oi? {eT subeotad o3 x0 yah py sd Ta 3 ise - * rats 
sn tam et a oat 

; ee Reno = EN ah0. = Lad 3 a x 


tat tac) ence - xix) Tree “ii si _ ie 
eorsixey bas. ores assem dsiw bedwetiitelh etiaannais wit >isosamy 
xb bejin ats 


-iieasi dtiw W Yo aclaulovaco iti) ety @) (3) 


II. SIMULATION 


The primary object of the simulation is to investigate the dis- 


tribution of the statistic B(n): 


[£ (x) - £(x) 1? 
B(n) -{ See G So? dx , 


over a suitable range of integration. We have broken the simulation 


into the following cases: 


A. UNIFORM (0,1) CASE 
Since Ihr O<x<l 
£(x) = 


oO, otherwise 


B(n) becomes 
1-b (n) 
B(n) -f [£(x)-117 ax . 
b(n) 
We integrate from b(n) to 1-b(n) instead of from 0 to 1, simply 
because we want to avoid the marked bias near 0 and 1. It can easily 


be shown that fF (x) is biased in the end regions, but not otherwise. 
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Thus the mean and variance of £ (x) in this case is: 


1 
ree Se x-y 
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F 7 
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This is only true for b(n) < x < l-b(n). Also 
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The integration was performed in the simulation using Simpson's 
rule, dividing the interval of integration into 100 equally spaced 
sub-intervals. The results were found to be satisfactory in the 
sense that the value of the integral for 100 sub-intervals differs 


little from those with more sub-intervals. 
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To evaluate £{*)» use the following formula: 


where the weight function W(u) is triangular, i.e.: 


1-|ul, if) [al <2 
W(u) = 
(0) A otherwise 


X50 j=l, ..., n are independent and identically distributed random 
variables. In this first case, they aré uniform (0,1) random variables. 
Also b(n) is the bandwidth. We employed three bandwidths in the 
uniform (0,1) case, namely: 3/vn, 1jva, l/n, one at a time. For each 
bandwidth, we used five different sample sizes. i.e. 100, 200, 500, 
1000, 1500, and for each sample size, we used five different random 
number seeds, i.e. simulated five times. The reason for this is to 
have the abiTity to assess the variability of the simulation results 
from the sections. The random number seeds are taken from the imme- 
diate previous "1500 sample size" output, except the first random 
number seed which is arbitrary. So, there are five sets of random 


numbers which are entirely different from one another. The following 


table should help to clarify the scheme: 
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TABLE I 


Bandwidth Sample size, n No. replications # Sections 


For each sample size, n, repeat five times using different random 
number seeds. 
For the uniform (0,1) case, there are 75 outputs altogether. In 
- each single simulation (one computer output), we evaluated 400 B(n)'s, 
say B5(n), j=l, .--, 400. The output includes the esinouines 
1. A bar histogram of the 400 £(n)'s. 
2. SNAP/IEDA graphs: 
a. A refined solid histogram. 
b. Graph of empirical log-survivor function, %n(1-j/(N+1)) vs 


Bi5) (me where B 


(n) is an ordered value. 
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c. Empirical CDF graph i/(n+l) vs B (5) (m). where Bis) ™ 
denotes the ordered observed B5(n)"s. 

d. Normal probability distribution with estimated mean and 
variance. Norm vs e.c.d.f. of Bj (n), where Norm is computed as 


follows: Let 


and s7(B(n)) = sample variance of the B;(n). Then 


7 Bm - i 
vetted Be [ S(B(n)) | 


e. Normal probability plot. Norm vs i/(N+l). 
A sample output is included in the "Computer Output" section 
(page 30 - 35). 

For each bandwidth and each sample size, there are five computer 
outputs, or 5 x 400 = 2000 B,(n)'s. With this 2000 sample size, we 
constructed a histogram again. Also we used the Rosenblatt smooth- 
ing technique to get a smooth ereiey. cunction estimate for B(n). 
After a number of trials, we found that the bandwidth needed in the 
Rosenblatt estimator has very much to do with the smoothing and how 
"good" the density looks. Fortunately it is robust to a certain 
extent. This graph can give us a better picture of the true distri- 
bution than the classical histogram. The last graph of this complete 


output for sample size 2000 is the empirical log survival function; 


£n(1 - i/(N+1)) vs 8B (n) e 


(3) 
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There are 15 "2000 sample size" computer outputs in the uniform (0,1) 

‘case. (Since each bandwidth has five "2000 sample size" outputs.) 
The complete computer outputs for each bandwidth and each sample 

size in the uniform (0,1) case are included in the "Computer Output" 


section (page 42 - 86). 


B. THE CAUCHY DISTRIBUTION 
The next case investigated is the Cauchy distribution, again with 


a triangular weight function and with the range of integration (-3,3): 


+3 2 
(£ (x) - £(x)] 
B(n) = = En A 
= 
The interval of integration comprises 80% of the probability mass for 
this distribution. The integration is carried out, in this case, 
using Simpson's rule with 600 sub-intervals. This grid was used after 
examining the 100, 300, 600, 900 sub-interval cases to ensure that a 
fine enough grid is used in the numerical integration. 
The Cauchy density function is 


£ (x) - ° 


T(1 + =) 


The density function estimator is the same as in the uniform case: 


n 
£ 0x) = ao dy wf} 


The weight function W(u) is again triangular. Thus the only differ- 


ence in this formula from the previous case is that Xie j=l, ices i 
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are i.i,d Cauchy random variables. Note also that for finite n the 
estimator will have a bias component, this bias component usually 
decreasing, for given n,as bandwidth decreases. In that case though 
the (pointwise) variance of £ fs) increases and possibly also the 
variance of 8(n). Note that &(n) and its distribution is a function 
of both sample size and bandwidth. 

The bandwidth b(n)'s used here are 1/vn, 3/vn, 20/Yn, the last 
one representing a case in which bias in the estimate £ Lf) plays a 
major role in ane distribution of 8(n). 

The rest of the simulation scheme is exactly the same as for the 
uniform (0,1) case. 

A sample output for the 400 B,(n) case is included in the "Com- 
puter Output" section (page 36 - 41). Also the complete computer 
outputs (2000 replications) for each bandwidth and each sample size 
in the Cauchy case are included in the "Computer Output" section 


(page 87 -131). 


C. HAND PLOT GRAPHS 

Besides the computer work, there are six hand plot graphs. These 
are the §(n) quantiles versus the standard normal quantiles, They 
are included in the Appendix, to illustrate the rate of convergence 


of the distribution of B(n) to its asymptotic normal distribution. 
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III. RESULTS 


The asymptotic result, conjectured by Rosenblatt, is that by 
using a triangular weight function when estimating the uniform 


density function on (0,1) 
1=b(n) 


v1/2 2 2 
{b(n) } nb (n) | £(x)-1|“dx-(1-2b(n)) J wW(u) “du 
b(n) 


is asymptotically normally distributed with mean zero and variance 


aw '4) (9) (1-2b(n)) 


2/9 


as n> if nb(n) > © and b(n) = o(n_ ). Here n is the sample 


size and b(n) the bandwidth. For the triangular weight function 


J W(u) 2du = 2 


=—co 


(4) 


and W (0) the 4th convolution of W with itself at zero, is 302/630. 


From the above expressions, we get: 


1-b(n) 
= -1|2 2 (1-2b(n)) 
E(B(n)) = -[ [ |£, (x) 1| «| ee ee 
b(n) 
1-b(n) ae 
var (B(n)) = var | £, 00-1 |7ax a cL MATCIIUEE ICAP 
n b(n) 
b(n) 


A comparison of the simulated values with that of the conjectured 
is tabulated below for the different cases investigated when f(x) is 


a uniform density: 
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TABLE II 


Comparison of estimated mean values and 
asymptotic mean values of B(n) for dif- 
ferent bandwidths and sample sizes. 


E(B (n))/(1-2b(n) ) 
b(n)=3/Yn | E(B(n)) Conjectured Computer output 


3000 -0089 20222 0127 
2121 - 0090 20157 0109 
21342 0073 0099 0075 
-0949 -0057 0070 -0058 
20775 0048 0057 -0051 


E(B(n))/(1-2b(n) ) 
b(n)=1/vn E(B(n) ) Conjectured Computer 


TABLE III 


Comparison of estimated standard deviation 
values and asymptotic standard deviation 
values of 8(n) for different bandwidths and 
sample sizes. 


0(B(n))/(1-2b(n) ) 
o(B(n) ) Conjectured | Computer output 


- 3000 0113 
2121 -0081 
~1342 -0046 
0949 -0029 
-0775 20022 


o(B(n))/(1-2b(n) ) 
b(n)=1/v7n o(B(n)) Conjectured | Computer output 
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Eacecialiy for small bandwidth the agreement between the asymptotic 
and simulated variances is very good even for small n(n=100). The 
same is true for the expected value, although convergence is slower 
than for 0 and again slower for large bandwidth. 

The simulation reveals that the distribution of the B(n) statistic 
converges to normality very slowly (cf. the six eae capa in the 
Appendix). For sample sizes 100, 200, 500, 1000, 1500, the histograms 
are all skewed to the right; the empirical log-survivor function plots 
are concave downward; and none of the normal probability plots is 
approximately straight lines (refer to "Computer Output" section). 

It is very suggestive that the B(n) statistic is approximately Gamma 


(8,k) distributed where the density of the Gamma distribution is: 


k . = 
eo Wick aie 
(k-1) ! 
and the mean and variance are: 
- 2 
E(X) = 0; Var(X) = kO"” . 


Therefore, for each underlying distribution, each bandwidth, and 
each sample size, we used the 2000 evaluated B(n)'s to estimate the 


6's and k's, and tabulate them as follows: 
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TABLE IV 


Estimated parameters for fitted Gamma distribution for 
B(n) in the case where the densities are uniform (0,1) 
and Cauchy, with a triangular weight function. 


(1) Bandwidth = 1/n Underlying distribution: Uniform (0,1) 


Sample sizes and estimated 9's and k's 


iy 100 200 1500 (std. dev.) 


S| 6520] 6638] .6612| 6621] 6645 = 0045 | 


k 39.48 38.46 32.95 30.86 31.06 + 2.91 


(2) Bandwidth = 1/vn Underlying distribution: Uniform (0,1) 


Sample sizes and estimated 8's and k's 


| | 100 200 500 1000 1500 (std. dev.) 
i 0565} .0420{ .0279| .0201 -0166 + .0003 
3.62 5.56 8.631 12.61 16.92 + 1.45 


Sample sizes and estimated 6's and k's 


Underlying distribution: Cauchy 


1500 (std. dev.) 
8.3921] .2800| .1764|_.1258] 1029 +0004 
k 21.02 31.49 59.34 79.19 100.83 + 4.37 

(3) Bandwidth = 3/vYn Underlying distribution: Uniform (0,1) 


sae sizes and estimated 9's and k's 
100 200 500 1000 } 1500 (std. dev.) 


Mies | daa 
1.59 2.43 3.65 5.03 + 0.37 


Sample sizes and estimated 6's and k's 


Underlying distribution: Cauchy 


1500 (std. dev.) 
a  —— 
k 8.62 12.30 20.32 29.14 33.98 1.55 

(4) Bandwidth = 20/vn Underlying distribution: Cauchy 


le sizes and estimated 9's and k's 
100 (std.dev.) | 200 500 1000| 1500 (std. dev.) 


}6|.o5se + .0008 |.0262| .o110|.0066].0052 + .0001 


Samp 


where 6 = mean > 
~ 2 
k = 1/(estimated coefficient of variation) 
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The standard deviation estimates for the parameter estimates are 
obtained from the five identical simulations using different random 
number seeds, each generating a sample size of 400, as detailed in 


Table I. 
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IV. CONCLUSIONS 


Although the conjecture proposed by Rosenblatt is that the dis- 
tribution of the statistic B(n) is asymptotically normal, this is 
certainly not the case for moderate sample sizes such as those used 
in the simulations. Based on this fact, the conjecture has little 
practical value, but the fact that all cases examined give rise to 
an approximate Gamma distribution gives some hope for stability in 
the results. 

However, the method proposed by Rosenblatt for estimating density 
functions is no doubt plausible as can be seen from the Rosenblatt 
type estimates used in this study. More effort in this area is 
needed, particularly to relate the Gamma distribution parameters to 
underlying factors. Since we realize that for moderate sample sizes, 
the statistic B(n) is not distribution free, and the bandwith b(n) 

(a function of sample size) plays a crucial role in estimating the 
density Fencevon and also the shape of the weight function W(u) is 
a major factor for evaluating £ fx) the future work may be to find 
the best combination of the above three factors by simulating the 
distribution of the statistic B(n) for different underlying distri- 
butions, bandwidths, sample sizes and weight functions and tabulate 
the resulting quantiles. Then we may be able to put the statistic 


B(n) to actual use. 
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APPENDIX 


Included in this Appendix are six plots of sample size versus 
estimated quantiles of the statistic B(n), normalized by subtract- 
ing the estimated mean and dividing by the estimated standard 
deviation. Also plotted as straight lines are the corresponding 
quantiles of the unit normal distribution. 

The 17 quaneties, extracted from each 2000-B(n) sample, are 
0.010, 0.020, 0.025, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 
0.80, 0.90, 0.95, 0.975, 0.980, 0.990. 


Other data for each plot are self-explanatory. 
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COMPUTER OUTPUT 


There are altogether 102 sheets of graphs. All of them are self- 


explanatory, except that the letter "X" in the graphs represents the 


statistic B(n). 


The following is a list of index of the graphs. 


36 - 41 


Description 


sample computer output for the case: 
underlying distribution-uniform (0,1); 
underlying distribution sample size n=1500; 
bandwidth 3/vn; 

statistic B(n) sample size m=400. 


sample computer output for the case: 
underlying distribution-Cauchy; 
underlying distribution sample size n=100; 
bandwidth 1/vn; 

statistic B(n) sample size m=400. 


The rest of the computer graphs in this section are all for statistic 


B(n) sample size m=2000. The page index is as follows. . 
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Underlying distribution 

distribution sample size Bandwidth 
u(0,1) 100 3/vn 
U(0,1) 200 . 3/vn 
u(0,1) 500 3/vn 
u(0,1) 1000 3/vn 
u(0,1) 1500 3/vn 
u(0,1) 100 1/vn 
u(0,1) 200 1/vn 
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